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An algorithm for calculating generalized fractal dimension of a time series using the general in- 
formation function is presented. The algorithm is based on a strings sort technique and requires 
0(N log 2 N) computations. A rough estimate for the number of points needed for the fractal di- 
mension calculation is given. The algorithm was tested on analytic example as well as well-known 
examples, such as, the Lorenz attractor, the Rossler attractor, the van der Pol oscillator, and the 
Mackey-Glass equation, and compared, successfully, with previous results published in the literature. 
The computation time for the algorithm suggested in this paper is much less then the computation 
time according to other methods. 
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I. BACKGROUND 



In the recent decades the study of chaos theory has gath- 
ered momentum. The complexity that can be found in 
many physical and biological systems has been analyzed 
by the tools of chaos theory. Characteristic properties, 
such as the Lyapunov exponent, Kolmogorov entropy and 
the fractal dimension (FD), have been measured in ex- 
perimental systems. It is fairly easy to calculate signs of 
chaos if the system can be represented by a set of non- 
linear ordinary differential equations. In many cases it is 
very difficult to built a mathematical model that can rep- 
resent sharply the experimental system. It is essential, 
for this purpose, to reconstruct a new phase space based 
on the information that one can produce from the system. 
A global value that is relatively simple to compute is the 
FD. The FD can give an indication of the dimensionality 
and complexity of the system. Since actual living bio- 
logical systems are not stable and the system complexity 
varies with time, one can distinguish between different 
states of the system by the FD. The FD can also deter- 
mine whether a particular system is more complex than 
other systems. However, since biological systems are very 
complex, it is better to use all the information and details 
the system can provide. In this paper we will present an 
algorithm for calculating FD based on the geometrical 
structure of the system. The method can provide impor- 
tant information, in addition, on the geometrical of the 
system (as reconstructed from a time series). 

The most common way of calculating FD is through the 
correlation function, C q (r) (eq. ^|). There is also other 
method of FD calculation based on Lyapunov exponents 
and Kaplan- Yorke conjecture Q (eq. 27). However, the 
computation of the Lyapunov exponent spectrum from a 
time series is very difficult and people usually try to avoid 



this method^. The algorithm which is presented in this 
paper is important, since it gives a comparative method 
for calculating FD according to the correlation function. 
The need for a additional method of FD calculation is 
critical in some types of time series, such as EEG series 
(which are produced by brain activity), since several dif- 
ferent estimations for FD were published in the literature 
A comparative algorithm can help to reach final 
conclusions about FD estimate for the signal. 



A very simple way to reconstruct a phase space from 
single time series was suggested by Takens 0. Giving a 
time series x j, i = 1 
phase space in the following way 



N p , we build a new n dimensional 



Vo = {x(t ),x(t + t), x(t Q + 2t), . 
yi = {x(ti),x(ti + T),x(h + 2t), 



,x{t + (n-l)r)} 
,ar(ti + (n-l)r)} 



y N = {x(t N ),x(t N + r),x(t N + 2t), . . . , x(t N + (n - l)r)} 



U — to + iAt t = mAt 

N = N p — (n — l)m m = integer, 



(1) 



where At is the sampling rate, r corresponds to the in- 
terval on the time series that creates the reconstructed 
phase space (it is usually chosen to be the first zero of the 
autocorrelation function Q , or the first minimum of mu- 
tual information ^ ; in this work we will use m instead of 
t as an index), and N is number of reconstructed vectors. 



"The basic difficulty is that for high FD there are some expo- 
nents which are very close to zero; one can easily add an extra 
unnecessary exponent that can increase the dimensionality by 
one; this difficulty is most dominant in Lyapunov exponents 
which have been calculated from a time series. 
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For ideal systems (an infinite number of points without 
external noise) any r can be chosen. Takens proved that 
such systems converge to the real dimension if n > 2D+1, 
where D is the FD of the real system. According to Ding 
et al [[l0| if n > D then the FD of the reconstructed phase 
space is equal to the actual FD. 

In this paper we will first (sec. [n|) present regular ana- 
lytic methods for calculating the FD of a system (gen- 
eralized correlation method and generalized information 
method). An efficient algorithm for calculating the FD 
from a time series based on st ring sorting will be de- 
scribed in the next section (sec. ID). The following step 
is to check and compare the general correlation methods 
with the general information method (see eq. (|^)) using 
known examples (sec. IV). Finally, we summarize in sec. 







II. 



GENERALIZED INFORMATION DIMENSION 
AND GENERALIZED CORRELATION 
DIMENSION 



A. Generalized information dimension 



D q = - lim 

e^o lne 



(4) 



for N — * co (N is the number of points). However, in 
practice, the requirement e —> is not achieved, and an 
average over a range of e is required. 





9 squares grid 

16 squeres grid 

























































































































The basic way to calculate an FD is with the Shannon 
entropy, ii(e)[[] [0, of the system. The entropy is just a 
particular case of the general information which is defined 
in the following way Jl2[ | : 

M(e) 



where we divide the phase space to M(e) hypercubes of 
edge e. The probability to find a point in the i th hyper 
cube is denoted by pi. The generalization, q, can be any 
real number; we usually use just an integer q. When we 
increase q, we give more weight to more populated boxes, 
and when we decrease q the less occupied boxes become 
dominant. For q = 1, by applying the l'Hospital rule, eq. 
(||) becomes : 

M(e) 

h[e) = - 2_, P< In Pi, (3) 

i=l 

where I\ (e) is referred to also as the Shannon entropy ]Tl| ] 
of the system. The definition of the general information 
dimension (GID) is : 



t Generally, one has to add a superscript, n, the embedding 
dimension in which the general information is calculated. At 
this stage we assume that the embedding dimension is known. 



FIG. 1. 2D Cantor set. Two different positions of the two 
dimensional grid, give completely different FD. 

In some cases this average is not sufficient because sev- 
eral values of I q {e) can be computed for the same e. To 
illustrate this we use a 2D Cantor set, presented in Fig. 
[l]. We start from a square. From 9 sub-squares we erase 
the 5 internal squares. We continue with this evolution 
for each remainder square. This procedure is continued, 
in principle, to infinity. The FD, Do, of the 2D Can- 
tor set is In 4/ In 3 (when q = 0, Iq is just the number 
of nonempty squares and Dq is the logarithmic ratio be- 
tween nonempty squares and e, where the square edge 
is normalized to one) as shown in Fig |l|. However, it is 
possible to locate the 2D grid in such a way that there 
are 16 nonempty squares, giving rise to Do = 2 In 4/ In 3, 
twice the time of real FD of the system. This illustra- 
tion shows that when e is not small, different positions 
of the grid can lead to different FD's. In this case it is 
clear that we have to locate the grid in such a way that 
minimum squares will be nonempty (it is easy to show 
that this claim for the minimum is true for every q). In 
general, we can say that one must locate the hyper-grid 
so that the general information is minimum : 

M(e) 

/,( e ) = min j— ln(X>?)> (5) 

5 i=l 

This proper location of the hyper-grid reduces the influ- 
ence of surface boxes that are partly contained in the 
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attractor. The requirement of a minimum gives a good 
estimate for the GID when e is not small. 



over all distances are associated with a loss in information 
based on the attractor structure. 



B. Generalized correlation dimension 



In 1983 Grassberger and Procaccia presented a new 
method for calculating Di (l3). According to this 
method, it is possible to calculate dimension just from 
the correlation between different points, without direct 
connection to the phase space, and therefore easy to use. 
Some years later, a more general correlation function was 
suggested by Pawelzik and Schuster O] : 



As we pointed out earlier, we usually have a limited num- 
ber of points, forcing us to calculate dimension at large r. 
Thus an error enters the calculation of FD. The minimum 
number of points needed for the FD calculation has been 
discussed in several papers and different constraints sug- 
gested (for example JItJ and (ll|). In this paper we will 
use the N m i n of the Eckmann and Ruelle [fL8| constraint 



Do < 



21ogioiV 
logio(i) 



(10) 



C q {r) 



N r N 

ly ly e 

i=l 1 7 = 1 



9-1 



(6) 



where N is number of points, point of the system 

and q is the generalization. <d(x) is the Heaviside step 
function 



Q(x) 



when x < 
when x > 0. 



(7) 



According to this method we have to calculate the gen- 
eralized probability to find any two points within a dis- 
tance r. This method is some kind of integration on eq. 
(|J). It is not necessary to compute the real distance (e.g. 
||x|| = \/ x\ + ■ ■ ■ + x 2 t where n is the phase space di- 
mension); it is equivalent to calculating the probability 
to find any two points in a hyper-box where one of them 
is located in the middle (e.g. ||x|| = maxK^,, \x.i\ [p^|). 
It is easier to compute the last possibility. For the special 
case of q = 1, cq. (|^) can be written (by applying the 
rHospital's rule 11) as : 



N r N 

lnCi(r) = — V In - Ve 
V ; N ^ N ^ 



Ft 



(8) 



under the following conditions 
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p = — < 1 



-N 2 p D » 1, 



(11) 



with reference to Grassberger and Procaccia method. 
Here, tq is the attractor diameter, and r is the maxi- 
mum distance that gives reliable results in eq. (Q) when 
q = 2. The normalized distance, p, must be small because 
of the misleading influence of the hyper-surface. A p too 
large (close to 1) can cause incorrect results since we take 
into account also hyper-boxes that are not well occupied 
because the major volume is outside the attractor. 

However, if one has a long time series, then it is not 
necessary to compute all iV 2 relations in eq. @. One 
can compute eq. (|^) for certain reference points such 
that the conditions in eq. ( pd] ) will still hold. Then, eq. 
ffl) can be written as follows : 



C q {r) 



N rcf 
Nref h 



Ndata -2W-1 



N r , 



E 01 



1 9-11 



T \Xi. s Xj\ 



(12) 



The generalized correlation dimension (GCD) has a sim- 
ilar definition to the GID (eq. (Q)) : 



D„ 



lim 

N— >oo;r— >0 



\nC q {r) 
lnr 



(9) 



Both GID and GCD methods should give identical re- 
sults. The GCD method is easy use and gives smooth 
curves. On the other hand, the method requires 0{N 2 ) 
computations^ Also, the smooth curves due to averaging 



*If one looks just at small r values, the method requires just 
O(N) computations Jl5[. 



where 



j\> w 



|^ data j 



N, 



ref 



For each reference point we calculate the correlation func- 
tion over all data points. The step for the time series is 
s. To neglect short time correlation one must also intro- 
duce a cutoff, w. Usually w ~ m (m corresponds to the 
r from (|]J)) p"9] ]. The number of distances \xi — Xj\ will 
be P = N re flN data - 2w - 1) instead of N 2 . The new 
form of eq. ([l0|) is: 



D < 



logio(^)' 



(13) 
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and the conditions ( |ll| ) become 



p=-«l 



2 P P » 1 



(14) 



Although we discussed here a minimum number of data 
points needed for the dimension calculation, one must 
choose p such that the influence from the surface is neg- 
ligible (the growth of the surface is exponential to the 
embedding dimension). That gives us an upper limit for 
p. On the other hand, in order to find the FD we must 
average over a range of r, giving us a lower limit for 
p; therefore N m i n is determined according to the lower 
limit, giving rise to larger amount of points. 



III. ALGORITHM 



A. Background 



They have proved that it is very efficient to divide the 
phase space into circles instead of squares in spite of the 
neglected area between the circles. However, this ap- 
proach is not suitable for higher phase space dimensions. 
The reason might be that the volume of the hyper-balls 
compared to the entire volume of the attractor decreases 
according to the attractor dimension, and in this way we 
lose most of the information that is included in the at- 
tractor, since the space between the hyper-balls is not 
taken into account. It is easy to show that the ratio be- 
tween the volume of the hyper-sphere with radius R, Vs, 
and the volume of the hyper-box with edge size of 2R, 
Vb (the sphere is in the box), is 



1 



Vs(2n) = 

V B (2n) 2-4-6--2n 



(£)"• 



for even phase space dimension, and, 



1 



Vs(2n - 1) _ 

V B (2n-l) 1 -3-5- --(211-1) V2 



(15) 



(16) 



In this section we describe an algorithm for GID 
method^]. The algorithm is based on a string sort and 
can be useful both for GID and GCD methods. 

In previous works there have been several suggestions to 
compute DIG |20| j25]|. The key idea of some of those 
methods p2| [24fl is to rescale the coordinates of each 
point and to express them in binary form. Then, the 
edge-box size can be initialized by the lowest value pos- 
sible and then double it each time. Those methods re- 
quire 0(N log N) operations. Another method |^l| uses 
a recursive algorithm to find the FD. The algorithm 
starts from a d dimensional box which contains all data 
points. Then, this box is divided into 2 d sub boxes, and 
so on. The empty boxes are not considered, and those 
which contain only one point are marked. This proce- 
dure requires only O(N) operations (for a reasonable se- 
ries length this fact does not make a significant difference 
f§)- 

As pointed out in ref. plj| , in spite of the efficiency of 
the above algorithm (speed and resolution), it is quite 
difficult to converge to the FD of a high dimensional sys- 
tem. We were aware to this difficulty, and suggested a 
smoothing term to solve it (as will be explained in this 
section). Basically, we allow any choice of edge length 
(in contrast to power of 2 edge size of the above meth- 
ods) and optimal location of the grid is searched. This 
procedure leads to convergence to the FD. 

One of the works that calculates GID was done by 
Caswell and Yorke pfj. They calculate FD of 2D maps. 



^Computer programs are available 

http:/ /faculty, biu. ac.il/ ~ ashkenaz/FDprog/ 



at 



for odd phase space dimension. It is clear that the ratios 
( |l5| ) and (|l^) tend to zero for large n. 



B. Algorithm 



Let us call the time series x(i). The form of the recon- 
structed vector, (|l|), depends on the jumping, to, that 
creates the phase space. We order the time series in the 
following way : 



x(0), x(m), x(2m), . . . , x((l 
x(l), x(l + to), x(l + 2m), . 



1)to), 

.,x(l + (l- l)m), 



x(m — 1), x(m — 1 + to), . . . , x{m — 1 + (I — 1)to), 



(17) 



where I = I — I . Let us denote the new series as Xj 
[i = . . . (lm — 1); we lost N p — Im data points). If 
we take n consecutive numbers in each row, we create 
a reconstructed vector in the n dimensional embedding 
dimension. 

At this stage, we fit a string (for each e) to the number 
series (0). One character is actually a number between 
. . . 255, and thus, it is possible to divide one edge of the 
hyper-box to 255 parts (we need one of the 256 characters 
for sign). The correspondence is applied in the following 
way. We search for the minimum value in $i series, and 
denote it as x m i n (similarly, we denote the maximum 
values as x max ). The corresponding character to Xi is : 



(18) 



4 



The value of e is in the range 

{5^ max %min ) 



255 



< e < (x„ 



(19) 



If we take now n consecutive characters (string) , we have 
the "address" of the hypercube in which the correspond- 
ing vector is found. 

Let us represent the string as follows : 
«o = yoyi ■ ..yn-i,si =yi...y n ,.. .,si- n = yi- n ■ ■ -yi-u 

Sl-n+l = Vl ■ ■ ■ yi+n-1, ■ ■ ■ , S2(/-n) + l = V2l-n ■ ■ ■ V21-1, 



s (m-l)(i-n+l) 



V (171-1)1 ■ ■ ■ J/(m-l)J+n-l) • ■ ■ i 

s m(Z-n+l)-l = Vml—n ■ ■ ■ Vml—1- 



(20) 



Obviously, we do not have to keep each string in the 
memory; we can keep just the pointers to the beginning 
of the strings in the characters series, yi. The number of 
vectors / strings is m(l — n + 1). 

As mentioned, each string is actually the address of a 
hyper-box (which the vector is in it) in a hyper-grid that 
cover the attractor. The first character is the position of 
the first coordinate of the hyper-box on the first edge of 
the hyper-grid. The second character denotes the loca- 
tion on the second edge, and so on. We actually grid the 
attractor with a hyper-grid so that any edge in this grid 
can be divided into 255 parts. Thus the maximal num- 
ber of boxes is 255", where n is the embedding dimension 
(there is no limit for n). Most of those boxes are empty, 
and we keep just the occupied boxes in the hyper-grid. 




sorting 



S 4 S,S 5 S„S () S,S ; ^s, 

s.b. s.b. 



s.b. = same box 



FIG. 2. Illustration of the GID algorithm. We take series 
containing 12 data points between zero and one. The param- 
eter values are: m = 3, e = 0.4 and n = 2. For simplicity we 
assume that the lowest character is 'a'. 



The next step is to check how many vectors fall in the 
same box, or, in other words, how many identical "ad- 
dresses" there are. For this we have to sort the vector 
that points to strings, Si, in increasing order (one string is 
less then other when the first character that is not equal 
is less then the parallel character; e.g., 'abcce'<'abcda' 
since the fourth character in the first string, c, is less then 
the fourth character in the second string, d). The above 
process is illustrated in Fig. |[ 

The most efficient way to sort N elements is the "quick 
sort" (2^^]. It requires just 0(N\og 2 N) computations. 
However, this sort algorithm is not suitable for our pro- 
pose, since after the vector is sorted, and we slightly in- 
crease the size of the edge of the hyper-box, e, there are 
just a few changes in the vector, and most of it remains 
sorted. Thus, one has to use another method which re- 
quires less computations for this kind of situation, be- 
cause the quick sort requires 0(N log 2 N) computations 
independently of the initial state. The sort that we used 
was a "shell sort" j2f|, which requires 0(N 3 / 2 ) compu- 
tations in worst case, around 0(N 5 / 4 ) computations for 
random series, and O(N) operations for an almost sorted 
vector. Thus, if we compute the information function for 
to different e values 0{N\og 2 N + mN) operations will 
be needed. 

After sorting we count how many identical vectors there 
are in each box. Suppose that we want to calculate GID 
for embedding dimensions n m .;„ . . . n maxi and generaliza- 
tions q-min ■ ■ ■ <lmax ■ We count in the sorted pointers vec- 
tor identical strings containing n ma x characters. Once 
we detect a difference we know that the coordinate of 
the box has changed. We detect the first different loca- 
tion. First difference in the last character of the string 
means that it is possible to observe the difference just in 
the n^ aa . embedding dimension, and in lower embedding 
dimension it is impossible to observe the difference. If the 
first mismatch is in the n th character in the string, then 
the box changes for embedding dimensions higher than 
or equal to n. We hold a counter vector for the different 
embedding dimensions, and we will assign zero for those 
embedding dimensions in which we observed a change. 
In this stage, we have to add to the results matrix the 
probability to fall in a box according to eqs. (j|), (|^). 
We continue with this method to the end of the pointer 
vector and then print the results matrix for the current 
e and continue to the next e. 

It is easy to show that one does not have to worry about 
the range of the generalization, q, because for large time 
series with 10 6 data points the range of computational 
q's is -100. .100. 

The method of string sorting can be used also for calcu- 
lating GCD. Our algorithm is actually a generalization of 
the method that was suggested by Grassberger ]l5|] . The 
algorithm is specially efficient for small r (r < j^D where 
D is the attractor diameter) and it requires 0(N log 2 N) 
computations instead of 0(N 2 ) computations (according 
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to the Grassbeger algorithm it requires O(N) computa- 
tions instead 0(N 2 ) computations). 

The basic idea of Grassberger was that for small r max 
values, where r max is the maximum distance for which 
the correlation function is computed, one does not have 
to consider distances larger then r max , distances which 
require most of the computation time. Thus, it is enough 
to consider just the neighboring boxes (we use the defini- 
tion distance ||a?|| = maxi<K„ \xi\ in eq. (||)), for which 
their edge is equal to r max , since distances greater then 
Tmax are not considered in eq. (g). If, for example, one 
wants to calculate the correlation of r < r maa 

' (IT 



r D 
then 



(in accordance with conditions ( |10D ( Jl l[ ) (|13|) (|14|)) 
(in 2D projection) it is enough to calculate distances in 
9 squares instead 100 squares (actually, it is enough to 
consider just 5 squares since eq. (||) is symmetric). Ac- 
cording to Grassberger, one has to build a matrix for 
which any cell in it points to the list of data points lo- 
cated in it. It is possible to generalize this method to a 
3-D projection matrix. 

The generalization to higher dimensions can be done very 
easily according to the string sort method. As a first step 
one has to prepare a string yi (^), which is the same 
operation as gridding the attractor by a hyper-grid of 
size r. The next step is to sort strings Si (|20|). For each 
reference point in eq. ( |l2| ) the boxes adjoining the box 
with the reference point in it should be found. Now, it is 
left to find the distances between the reference point to 
other points in neighboring boxes (we keep the pointer 
vector from string yi to the initial series and vice versa), 
and calculate the correlation according to eq. (|l2|). 

The algorithm described above is especially efficient for 
very complex systems with high FD. For example, for 
EEG series produce by brain activity, it well known (ex- 
cept for very special cases such as epilepsy jlj|) that the 
FD is, at least, four. In this case n m i n — 4, and instead 
of calculating distances in l A boxes (I is the ratio between 
the attractor diameter and r) it is sufficient to calculate 
distances in 3 4 = 81 boxes. If, for example, / = 9, just 
3 4 /9 4 = gy distances must be computed (or even less if 
one take into account also the symmetry of eq. (|l2|)). In 
this way, despite the 0(Nlog 2 N) operations needed for 
the initial sort, one can reduce significantly the compu- 
tation time. 



C. Testing of GID algorithm and number of data 
points needed for calculating GID 

Let us test the algorithm of GID on random numbers. We 
create a random series between to 1 (uniform distribu- 
tion) . We then reconstruct the phase space according to 
Takens theory (Q). For any embedding dimension that 
we would like to reconstruct, the phase space we will get 
a dimension that is the same as the embedding dimen- 
sion, since the new reconstructed vectors are random in 



the new phase space, and hence fill the space (in our case, 
a hyper-cube of edge 1). 

For embedding dimension 1 (a simple line), if the edge 
length is e, then the probability to fall in any edge is also 
e. The number of edges e that covers the series range 
[0,1) is L i J . The probability to fall in the last edge is 
l%e (the residue of the division). For generalization q 
one gets : 



M(e) 



p?= [l/e\e ci + {l%e) q . 



(21) 



For embedding dimension 2, one has to square (|21|) since 
the probability distribution on different sides is equal. 
For embedding dimension n eq. (||) becomes : 



1-9 



ln[|l/eje 9 + (l%e) 



gin 



(22) 



Opening eq. ( p2| ) according to the binomial formula will 
give all different combinations for the partly contained 
boxes. Notice that this grid location fulfills requirement 
(o) for the minimum information function. 



-Ue) 




-10.0 ' 



-15.0 



-5.0 



-4.0 



-3.0 -2.0 

Ine 



-1.0 0.0 



FIG. 3. Dimension calculation for random series by GID 
algorithm and by the theory. 

In Fig. U we present —/2(e), both, according to the GID 
algorithm, and according to eq. (22). There is full cor- 
respondence between the curves, and, as we find, it is 
possible to see that the slope (dimension) in the embed- 
ding dimension n is equal to the dimension itself. When 
—/2(e) ~ —10, the curves separate. There is lower con- 
vergence (when —/2(e) ~ —12.5) since the edge is too 
small, and the number of reconstructed vectors is equal 
to the number of nonempty boxes. 
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The "saw tooth" that seen in Fig. [| caused by partial 
boxes that fall on the edge of the hyper-box. The "saw 
tooth" appears when there is an integer number of edges 
that covers the big edge. The slope in the beginning of 
every saw tooth is equal to zero (that comes from the 
extremum condition of (H^)). The "saw teeth" become 
smaller when e decrease, since the relative number of 
boxes that fall on the surface in small compared to the 
entire number of boxes. 

It is possible to calculate roughly the number of points 
needed for the dimension calculation in a homogeneous 
system. The separation point of Fig. |^ is approxi- 
mately around —I^ie) w —10, for every embedding di- 
mension. Thus, the number of hyper-boxes at this point 
is e 10 (since —10 = In ( y M(e)p 2 ) where p ~ jj^ in 
our example). The graph come to a saturation around 
— 7 2 (e) ~ —12.5, giving rise to e 12 5 data points. The 
number of points per box is just, e 2 5 w 12. Thus, gener- 
ally, when there are less then 12 points per box, the value 
of I q (e) is unreliable. Let us denote the minimum edges 
needed for the calculation as m m i n (m m i n is just |l/ e J> 
and thus one can define a lower value for e) . The number 
of hyper-boxes in embedding dimension n is Thus, 
the minimum number of points needed for computing the 
FD of an attractor with an FD D, is: 

N = 12m£ to . (23) 

This estimate is less then what was required by Smith 
p7[ (A2 D ) and larger then the requirement of Eckmann 
and Ruelle jl8) (if we take, for example, m m j n — | = 

10 then according to ( |l0|) 10 D / 2 points is needed and 
according to ( p3| ) 12 x 10 points is needed). However, 
as we pointed out earlier, ( p3| ) is a rough estimation, and 
one can converge to a desired FD even with less points. 



IV. EXAMPLES 

In this section we will present the results of both GID and 
GCD on well known systems, such as the Lorenz attrac- 
tor and the Rossler attractor. The FD of those results 
are well known by various methods, and we will present 
almost identical results achieved by our methods. In the 
following examples, we choose one of the coordinates to 
be the time series; we normalize the time series to be in 
the range of ... 1. 



A. Lorenz Attractor 

The Lorenz attractor p§] is represented by a set of three 
first order non-linear differential equations : 



§ = o-(y-*) 

— = —xz + rx — y (24) 
at 

dz 

— = xy — bz. 
dt y 

Although it is much simpler to find the FD from the real 
phase-space (x, y, z), we prefer the difficult case of finding 
the FD just from one coordinate to show the efficiency of 
Takens theorem. We choose the z th component to be the 
time series from which we reconstruct the phase space 
(the parameter values are: a = 16, r = 45.92, 6 = 4). 
We took 32768 data points and m — 6 (the jumping 
on the time series which creates the reconstructed phase 
space; the time step is At = 0.02). The generalization 
that we used was q = 2. The embedding dimensions were 
n = 1..7. The FD of Lorenz attractor is 2.07 ( |§9) and 
others). 

In Fig. U we show the results of GCD and GID methods. 
As expected from the GCD curves, in Fig. ||]a we see very 
smooth curves, for which the slopes (in the center parts 
of the curves) converge approximately to dimension 2.08. 
In Fig. we see some jumping in the GID curves. It 
can be clearly seen that for small e, one can not see the 
jumping while for large e the non-monotonicity is much 
stronger. This is typical behavior for GID curves; large e 
(or larger box edge) reflects more variability of —/2(e) to 
the hype- grid location. To smooth the curves in order to 
get reliable results, it is necessary to find a proper loca- 
tion of the hyper-grid that gives minimum general infor- 
mation I%{e) (or maximum —12(e)). We have performed 
7 comparisons to find that proper location. As for GCD, 
we also fitted approximate linear curves (dashed lines) 
which lead to dimension 2.09 (Fig. ^c), and thus, agree 
well with previous results. To make sure of the validity 
of the dimension results we add here (Fig. |^d) graphs of 
successive slopes of curves in Fig. ^c, and conclude with 
dimension results versus embedding dimension according 
to several generalizations (q — 2. .5, Fig. ||e). In both 
graphs, the approximate dimension is 2.03. 
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FIG. 4. Correlation and information curves of recon- 
structed Lorenz attractor. a. Correlation curves; the esti- 
mate slopes are added, b. Unsmoothed information curves, 
c. Smoothed information curves; the estimated linear curves 
(dashed line) and their slope is added, d. The slopes of —/2(e) 
(Fig. |^c). The dashed line represents the estimated FD. e. 
FD of Lorenz attractor for different embedding dimensions, 
n, and for different generalizations, q. 



B. Rossler Attractor 



Let us examine briefly our second example, the Rossler 
attractor pOfl . The differential equations which generate 
the trajectory are : 



dx 

~dt 
dy 

dt 

dz 

~dt 



-z-y 
x + ay 
b + z(x — c). 



(25) 
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FIG. 5. a. Unsmoothed graphs of Rossler attractor. b. 
Smoothed information graphs; the FD estimation is added. 



The parameter values are : a = 0.15, b = 0.2, c = 10. 
The time step is t = 0.2 and m = 8. We use the y direc- 
tion to be the time series. The number of data points is 
N = 16384. In Fig. | we show the results of GID method 
of embedding dimension n = 1..7. There are many more 
discontinuities in the unsmoothed curves (Fig. ||a) com- 
pared to the Lorenz attractor, reflecting a sensitivity to 
the grid location. The smooth curves in Fig. |5|b are pro- 
duced after 20 comparisons. Again, we approximate the 
slopes of central linear part of the curves in two differ- 
ent ways and find dimension 2.05, which is in agreement 
with dimension 2.03 calculated by the use of the GCD 
method, and through the Lyapunov exponent p9|. 



C. Van der Pol Oscillator 

Another example that will be tested is the van der Pol 
oscillator |5l] . This system was investigated in the begin- 
ning of the century in order to find a model for the heart 
beat variability (among other uses of this equation). The 
behavior of the system (with parameter values that will 
be used) is chaotic although it looks quite periodic. Thus, 
we expect a smooth behavior of the geometrical structure 
of the attractor. The equation of motion is : 

x — a(l — x 2 )x + kx = f cos fit. (26) 
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The parameter values are : O = 2.446, a = 5, k = 1 
and / = 1. Habib and Ryne |52j and others j33| have 
found that the Lyapunov exponents of the system are 
Ai « 0.098, A 2 = and A 3 « —6.84, and thus, according 
to the Kaplan and Yorke formula 



i+il 



0.098 
6.84 



2.014 



(27) 



where j is defined by the condition that A, > and 



Si=i < 0- The fact that the system has a "periodic" 
nature is reflected in this very low FD (as known, the 
minimal FD of a chaotic system is 2). 
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dimensional dynamics. A delay differential equation of 
the type 



dx(t) 
dt 



F(x(t),x(t-T)) 



(28) 



belongs to this class. In order to solve this kind of equa- 
tion, one needs to initiate x(t) in the range < t < r, and 
then, by assuming that the system can be represented by 
a finite set of discrete values of x(t), it is possible to eval- 
uate the system dynamics step by step. The solution is 
considered to be accurate if one converges to the same 
global properties, such as, Lyapunov exponents and FD, 
by different methods p4|], and without dependence on 
the number of intervals to which the function is divided 
p3[p|. Delay equations, such as eq. (p8|), describe sys- 
tems in which a stimulus has a delay response. There are 
many practical examples from control theory, economics, 
population biology, and other fields. 

One of the known examples is the model of blood cell pro- 
duction in pat ients with leukemia, formulated by Mackey 
and Glass @ : 



x(t) = r— -rz bx(t). 

W 1 + [x(t - t)] c V ; 



(29) 



Following refs. Q,^3[, the parameter values are : a — 
0.1, b — 0.2, and c = 10. We confine ourselves to r = 
30.0. As in ref. we choose the time series to be 

{x(t),x(t + r),x(t + 2r), . . as well as the integration 
method which is described in this reference. The length 
of the time series is N = 131072. 



FIG. 6. van der Pol oscillator : The smoothed information 
graphs, — Ji(e); The estimated FD is added (dashed lines with 
their corresponding slope). 

In Fig. ^ we present the calculation of D\ which we 
suppose to be close to Dl- The approximate slope lead 
to an FD of D\ ~ 2.02, which is in good agreement to 
Dl. Notice that minimization of I\{e) did not succeed 
for large e values although we used 7 grid comparisons; 
that fact can cause a small error in the slope estimation, 
since usually the slope is determined from the central 
part of the curves. 



D. Mackey-Glass Equation 

Up to now we examined the relation between GID and 
GCD on systems with finite dimension, for which their 
FD can be calculated by using methods based on the 
true integrated vector of the system. A class of systems 
which is more close to reality are systems with infinite 



"Note that the dynamics which is calculated by the use of 
different methods, or, equivalently, by a different number of 
discrete integration intervals, does not have to behave identi- 
cally. In fact, one can expect similar dynamics if the system 
is not chaotic, but, if the system is chaotic, it is sensitive 
to initial condition and to the integration method, and even 
moreover, it is sensitive to integration accuracy pB| , p6| . How- 
ever, the global properties of a system converge to the same 
values, by using different integration methods and using dif- 
ferent integration accuracy [M. 

t1 4n fact, the common procedure of reconstructing a phase 
space from a time series is described in eq. (|l|). According 
to this method, one has to build the reconstructed vectors by 
the use of a jumping choice m which is determined by the first 
zero of the autocorrelation function, or the first minimum of 
the mutual information function 0. However, we used the 
same series as in ref. [F3 in order to compare results. 
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FIG. 7. Mackey-Glass equation: a. The GCD graphs. 
One notices two regions of parallel curves, which lead to dif- 
ferent FD's. b. The successive slopes of a. c. The GID 
graphs, d. The successive slopes of c. 



The FD (D2) calculation of the Mackey-Glass equation 
is presented in Fig. 7L The embedding dimensions are, 
n = 1 ... 9. In Fig. [7a, the correlation function, C^r) is 
shown. One notices that there are two regions in each one 
of which there is convergence to a certain slope. These 
regions are separate by a dashed line. In Fig. [7]b, the 
average local slope of Fig. 0a is shown. One can identify 
easily two convergences, the first in the neighborhood of 
~ 3.0, and the second around ~ 2.4. Thus, there are two 
approximations for the FD, D2, pointing to two different 
scales. The first approximation is similar to the FD that 
was calculated in ref. Jl3]| . However, the GID graphs 
which are presented in Fig. 0c and d lead to an FD, 
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£>2 ~ 2.45, which seems to be very close to the second 
convergence of Fig. |^b. Notice that the convergence to 
D2 ~ 2.4 appears, in both methods, in the neighborhood 
of the same box size (~ 2.1). 



V. SUMMARY 

In this work we develop a new algorithm for calculation 
of a general information, I q (n), which is based on strings 
sorting (the method of string sort can be used also to 
calculate the conventional GCD method). According to 
our algorithm, one can divide the phase space into 255 
parts in each hyper-box edge. The algorithm requires 
0(N log 2 N) computations, where N is the number of 
reconstructed vectors. A rough estimate for the number 
of points needed for the FD calculation was given. The 
algorithm, which can be used in a regular system with 
known equations of motion, was tested on a reconstructed 
phase space (which was built according to Takens theo- 
rem). The general information graphs have non mono- 
tonic curves, which can be smoothed by the requirement 
for minimum general information. We examine our algo- 
rithm on some well known examples, such as, the Lorenz 
attractor, the Rossler attractor, the van der Pol oscillator 
and others, and show that the FD that was computed by 
the GID method is almost identical to the well known 
FD's of those systems. 

In practice, the computation time of an FD using the GID 
method, was much less than for the GCD method. For 
a typical time series with 32768 data points, the compu- 
tation time needed for the GCD method was about nine 
times greater than the computation time of GID method 
(when we do not restrict ourselves to small r values and 
we compute all N 2 relations). Thus, in addition to the 
fact that the algorithm developed in this paper enables 
the use of comparative methods (which is crucial in some 
cases), the algorithm is generally faster. 

The author wishes to thank J. Levitan, M. Lewkowicz 
and L.P. Horwitz for very useful discussions. 
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